Practical Session 9: The relation between two time series

Expected learning outcomes

By the end of this session, participants should be able to:

  • assess and interpret associations with external variables

  • identify and interpret effect modification between external variables and the outcome variable in surveillance data

Task 9.1

Using the aragon dataset, assess the effect of ambient temperature (using average weekly maximum temperature) on mortality in the autonomous community of Aragón. Is this effect uniform throughout the year?

Task 9.2 (Optional)

It has been argued by experts that low temperatures in winter might be “slow killers”; that is, very low temperatures in winter do not result in a peak in mortality on the same day/week as they are observed, but rather after some time has elapsed. On the other hand, very high temperatures in summer are “fast killers”; they are associated with peaks in mortality very fast, i.e. heat waves kill people fast. Can you find evidence for this for Aragón?

Help for Task 9.1

Open the dataset aragon.csv. There you will find mean maximum temperature and mortality data for the autonomous community of Aragon by week.

Show the code
aragon <- import(here("data", "aragon.csv")) 

Convert the data to a time series and plot both variables (weekly mean maximum temperature and mortality data).

Show the code
# Create ts dataframe with yearweek index
aragonz <- aragon %>% 
  mutate(
    year_week = make_yearweek(year = year, week = week),
    index = seq.int(from = 1, to = nrow(aragon))
  ) %>% 
  as_tsibble(index = index)

# Temperature plot   
aragonz_tmax_plot <- ggplot(data = aragonz) +
  geom_point(mapping = aes(x = year_week, y = tmax), colour = "blue", alpha = 0.5) +
  scale_y_continuous(limits = c(0, NA)) +
  scale_x_yearweek(date_labels = "%Y-%W", 
                   date_breaks = "1 year") +
  labs(x = "Week", 
       y = "Maximum Temperature") +
  tsa_theme

# Mortality plot
aragonz_deaths_plot <- ggplot(data = aragonz) +
  geom_point(mapping = aes(x = year_week, y = cases), colour = "green", alpha = 0.5) +
  scale_y_continuous(limits = c(0, NA)) +
  scale_x_yearweek(date_labels = "%Y-%W", 
                   date_breaks = "1 year") +
  labs(x = "Week", 
       y = "Weekly deaths counts") +
  tsa_theme

# Using patchwork to join the two plots
# The "/" determines an above/below display
aragonz_two_plot <- (aragonz_tmax_plot / aragonz_deaths_plot) 
aragonz_two_plot

Generate variables for sine and cosine for annual oscillation.

Show the code
aragonz <- aragonz %>% 
  mutate(sin52 = sin(2 * pi * date / 52),
         cos52 = cos(2 * pi * date / 52))

Fit a poisson regression model with a simple trend to the weekly number of deaths, accounting for seasonality.

Show the code
aragmodel1 <- glm(cases ~ index + sin52 + cos52,
                           family = "poisson",
                           data = aragonz)

aragmodel1 %>% 
  gtsummary::tbl_regression(
    exponentiate = T,
    estimate_fun = ~style_number(.x, digits = 4)
  )
Characteristic IRR 95% CI p-value
index 1.0001 1.0001, 1.0002 <0.001
sin52 1.0454 1.0364, 1.0545 <0.001
cos52 1.1094 1.0998, 1.1190 <0.001
Abbreviations: CI = Confidence Interval, IRR = Incidence Rate Ratio

Plot the predicted weekly number of deaths.

Show the code
ggplot(data = aragonz) +
  geom_point(mapping = aes(x = year_week, y = cases), alpha = 0.2) +
  geom_line(mapping = aes(x = year_week, y = aragmodel1$fitted.values), 
            colour = "darkorange", 
            alpha = 0.7, 
            lwd = 1.5) +
  scale_y_continuous(limits = c(0, NA)) +
  scale_x_yearweek(date_labels = "%Y-%W", 
                   date_breaks = "1 year") +
  labs(x = "Week", 
       y = "Weekly deaths counts", 
       title = "Deaths with fitted trend and seasonality") +
  tsa_theme

Discuss how your model fits the data.

Note: When you plotted the weekly number of deaths in Aragon, perhaps you noticed that mortality peaks in winter; however, there are smaller peaks during the summer, too. You can account for that in two different ways: a) include two sine/cosine functions in the model (one for 52-week cycles and one for 26-week cycles), or b) take into account the fact that temperature might be behaving differently as a risk factor in winter than it does in summer (sign of effect modification?).

You decide to perform your analysis twice; once for winter and once for the rest of the time. You define winter as the time between weeks 49 and 8.

Show the code
aragonz <- aragonz %>% 
    mutate(winter = (week >= 49 | week <= 8))

First run the model including winter as a main effect, and then including an interaction term with temperature.

Show the code
# Model with only main effect
aragmodel2 <- glm(cases ~ index + sin52 + cos52 + winter,
                           family = "poisson",
                           data = aragonz)

# Model with temperature interaction
aragmodel3 <- glm(cases ~ index + sin52 + cos52 + winter * tmax,
                           family = "poisson",
                           data = aragonz)

We compare both model’s AIC below.

Show the code
# Easiest, most direct way
AIC(aragmodel2); AIC(aragmodel3)
Show the code
# Fancier way using the library gtsummary
tbl_m2 <- aragmodel2 %>% 
  tbl_regression(exponentiate = TRUE,
                 estimate_fun = ~style_number(.x, digits = 4)) %>%
  add_glance_table(include = c(AIC, BIC, deviance))

tbl_m3 <- aragmodel3 %>% 
  tbl_regression(exponentiate = TRUE,
                 estimate_fun = ~style_number(.x, digits = 4)) %>%
  add_glance_table(include = c(AIC, BIC, deviance))

tbl_m2m3 <- 
  tbl_merge(
    tbls = list(tbl_m2, tbl_m3),
    tab_spanner = c("aragmodel2", "aragmodel3")
  )
tbl_m2m3
Characteristic
aragmodel2
aragmodel3
IRR 95% CI p-value IRR 95% CI p-value
index 1.0001 1.0001, 1.0002 <0.001 1.0001 1.0001, 1.0002 <0.001
sin52 1.0403 1.0311, 1.0495 <0.001 1.0621 1.0508, 1.0735 <0.001
cos52 1.0783 1.0651, 1.0917 <0.001 1.1721 1.1385, 1.2066 <0.001
winter





    FALSE

    TRUE 1.0673 1.0461, 1.0889 <0.001 1.2688 1.2020, 1.3392 <0.001
AIC 4,718

4,666

BIC 4,739

4,696

Deviance 1,006

950

tmax


1.0078 1.0054, 1.0103 <0.001
winter * tmax





    TRUE * tmax


0.9849 0.9806, 0.9893 <0.001
Abbreviations: CI = Confidence Interval, IRR = Incidence Rate Ratio

Note that winter * tmax tells R to include terms in the model for winter, tmax and for their interaction. If we just wanted to include the interaction term without its corresponding “main effects”, i.e. without winter or tmax, we would use winter:tmax instead.

Discuss the output of the two models. Does the interpretation change when winter is taken into account as an interaction term along with temperature?

Help for Task 9.2

You decide to check whether temperature has an effect on mortality with some lag. stats::lag is the default base R function to compute lags from the stats package. In the example below, we prefer to use one of the tidyverse package alternatives,namely dplyr::lag. We specify the package name here to avoid possible conflicts with other packages which have a lag function.

We are using the lag function from the dplyr package. For a given week, the value of tmax is lagged by 1 in respect to the previous week’s value.

Show the code
# Create lag variable
aragonz <- aragonz %>% 
  mutate(L1.tmax = dplyr::lag(x = tmax, n = 1))

Create a new model named aragmodel4 including the lag variable, and compare it with the previous model aragmodel3. What can you conclude from it?

Show the code
aragmodel4 <- glm(cases ~ index + sin52 + cos52 + winter * tmax + L1.tmax,
                           family = "poisson",
                           data = aragonz)

summary(aragmodel4)

Call:
glm(formula = cases ~ index + sin52 + cos52 + winter * tmax + 
    L1.tmax, family = "poisson", data = aragonz)

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)      5.087e+00  3.198e-02 159.076  < 2e-16 ***
index            1.246e-04  2.097e-05   5.942 2.81e-09 ***
sin52            5.679e-02  6.343e-03   8.954  < 2e-16 ***
cos52            1.504e-01  1.663e-02   9.043  < 2e-16 ***
winterTRUE       2.268e-01  2.795e-02   8.113 4.94e-16 ***
tmax             8.194e-03  1.304e-03   6.283 3.32e-10 ***
L1.tmax         -1.226e-03  1.146e-03  -1.070    0.285    
winterTRUE:tmax -1.433e-02  2.278e-03  -6.291 3.15e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 1708.45  on 518  degrees of freedom
Residual deviance:  943.72  on 511  degrees of freedom
  (1 observation deleted due to missingness)
AIC: 4654.1

Number of Fisher Scoring iterations: 4
Show the code
AIC(aragmodel3); AIC(aragmodel4)
[1] 4665.836
[1] 4654.111

Instead of running the model with an interaction term, you could run the analyses for winter and the rest of your dataset separately.

The interaction term was not significant, but you are still curious about the differential effect of temperatures on mortality in winter and the other seasons. You decide to split the data and run two different models for winter and non-winter mortality and temperatures, and compare them.

Think carefully about how you will compare them: are the AIC criteria still valid for this task?

Show the code
## winter subset
aragonz.winter <- 
    aragonz %>% 
    filter(winter == TRUE)

## not winter subset
aragonz.notwinter <- 
    aragonz %>% 
    filter(winter == FALSE)
Show the code
# Winter model
aragmodel5 <- glm(cases ~ index + sin52 + cos52 + L1.tmax,
                           family = "poisson",
                           data = aragonz.winter)

# Not winter model
aragmodel6 <- glm(cases ~ index + sin52 + cos52 + L1.tmax,
                           family = "poisson",
                           data = aragonz.notwinter)

# Table summaries
tbl_m5 <- aragmodel5 %>% 
  tbl_regression(exponentiate = TRUE,
                 estimate_fun = ~style_number(.x, digits = 4))

tbl_m6 <- aragmodel6 %>% 
  tbl_regression(exponentiate = TRUE,
                 estimate_fun = ~style_number(.x, digits = 4))

# Comparison table
tbl_m5m6 <- 
  tbl_merge(
    tbls = list(tbl_m5, tbl_m6),
    tab_spanner = c("Winter model", "Not winter model")
  )
tbl_m5m6
Characteristic
Winter model
Not winter model
IRR 95% CI p-value IRR 95% CI p-value
index 1.0002 1.0001, 1.0003 <0.001 1.0001 1.0001, 1.0002 <0.001
sin52 1.1944 1.1428, 1.2486 <0.001 1.0485 1.0348, 1.0624 <0.001
cos52 1.7229 1.4361, 2.0675 <0.001 1.1094 1.0785, 1.1412 <0.001
L1.tmax 0.9923 0.9879, 0.9967 <0.001 1.0031 1.0006, 1.0055 0.015
Abbreviations: CI = Confidence Interval, IRR = Incidence Rate Ratio

Discuss the output of the different models. Which one would you go for?